Solving Partial Differential Equations using Finite Elements

Oliver Rübenkönig

Wolfram Research

Why finite elements

Arbitrary shaped domains

Large class of PDEs

The coefficient form

Consider the single linear partial differential equation

m((∂^2)/(∂(t^2)))u+d(∂/(∂t))u+∇·(-c ∇u -α u + γ ) +β·∇u + a u = f on  Ω

Specific PDEs

Coefficient Form: m((∂^2)/(∂(t^2)))u+d(∂/(∂t))u+∇·(-c ∇u -α u + γ ) +β·∇u + a u = f

Laplace:                                                                       ∇·(-c ∇u -α u + γ) +β·∇u + a u = 0

Poisson :                                                                     ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Helmholtz:                                                                 ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Convection/Diffusion/Reaction:      ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Heat:                                        ((∂^2)/(∂(t^2)))u+(∂/(∂t))u+ ∇·(-c ∇u -α u + γ) +β·∇u + a u = f

Wave:                                     ((∂^2)/(∂(t^2)))u+(∂/(∂t))u+ ∇·(-c ∇u -α u + γ) +β·∇u + a u = 0

Heat conduction

∇·(-c ∇u)= f

PDE

c= 0.03;
pde=Div[-c Grad[u[x,y],{x,y}],{x,y}]==0

Region

r0=1/2;r1=1/10;d=3/20;
Ω=RegionDifference[RegionDifference[Disk[{0,0},r0],Disk[{-d,0},r1]],Disk[{d,0},r1]];

RegionPlot[True,{x,y}∈Ω]

Boundary conditions

(Subscript[Γ, D])={
DirichletCondition[u[x,y]==273,x^2+y^2==r0^2],
DirichletCondition[u[x,y]==323,(x+d)^2+y^2==r1^2],
DirichletCondition[u[x,y]==353,(x-d)^2+y^2==r1^2]
};

Solve

uif=NDSolveValue[{pde,(Subscript[Γ, D])},u,{x,y}∈Ω]

Plot

Plot3D[uif[x,y],{x,y}∈Ω]

cp=ContourPlot[uif[x,y],{x,y}∈Ω]//Quiet

Post process

{uifx,uify}=-c Grad[uif[x,y],{x,y}]

Show[
cp,
VectorPlot[{uifx,uify},{x,-r0,r0},{y,-r0,r0},VectorScale->Large,VectorStyle->Red,RegionFunction->Function[{x,y,z},Evaluate[RegionMember[Ω][{x,y}]]]]//Quiet]

Heat conduction with multiple regions

Coefficients can depend on space and time but not on the dependent variable.

∇·(-c(x,y) ∇u)= f

PDE

c=If[y≤sh,{{11.7,0.},{0.,11.7}},{{1.,0},{0.,1.}}];

ρc =10^-8.;ϵ0=8.86*^-12;
pde=Inactive[Div][-c.Inactive[Grad][u[x,y],{x,y}],{x,y}]==ρc/ϵ0;

pde=Diffusion[-c,u[x,y],{x,y}]==ρc/ϵ0

Region

sh=0.2;sh2=0.02;sw=0.3;
bmesh=ToBoundaryMesh["Coordinates"->{{0,0},{1,0},{1,sh},{1,1},{0,1},{0,sh+sh2},{sw,sh+sh2},{sw,sh},{0,sh}},"BoundaryElements"->{LineElement[{{1,2},{2,3},{3,4},{4,5},{5,6},{6,7},{7,8},{8,9},{9,1},{3,8}}]}];
bmw=bmesh["Wireframe"]

mesh=ToElementMesh[bmesh,"MaxCellMeasure"->0.001];
mesh["Wireframe"]

Boundary conditions

(Subscript[Γ, D])={
DirichletCondition[u[x,y]==0,x==1||y==1||y==0],
DirichletCondition[u[x,y]==10^3,0<x≤sw&&sh≤ y≤sh+sh2]
};

Solve

uif=NDSolveValue[{pde,(Subscript[Γ, D])},u,{x,y}∈mesh]

Plot

Show[ContourPlot[uif[x,y],{x,y}∈mesh],bmw]

Systems of equations

Consider a system of 2 equations for dependent variables u and v.

∇·(-(Subscript[c, 11]) ∇u -(Subscript[c, 12])∇v -(Subscript[α, 11]) u +(Subscript[α, 12]) v+ (Subscript[γ, 1]) ) +(Subscript[β, 11])·∇u +(Subscript[β, 12])·∇v + (Subscript[a, 11]) u+ (Subscript[a, 12]) v == (Subscript[f, 1])
∇·(-(Subscript[c, 21]) ∇u -(Subscript[c, 22])∇v -(Subscript[α, 21]) u +(Subscript[α, 22]) v+ (Subscript[γ, 2]) ) +(Subscript[β, 21])·∇u +(Subscript[β, 22])·∇v + (Subscript[a, 21]) u+ (Subscript[a, 22]) v== (Subscript[f, 2])

For a region Ω in (ℝ^n), each (Subscript[c, ij]) is an n×n matrix,and each of  (Subscript[α, ij]), (Subscript[β, ij]), and (Subscript[γ, i]) are vectors.

The Dirichlet boundary conditions are:

(Subscript[h, 11])u+(Subscript[h, 21])v=(Subscript[r, 1])
(Subscript[h, 21])u+(Subscript[h, 22])v=(Subscript[r, 2])

In this case, the generalized Neumann conditions are

(OverVector[n])·((Subscript[c, 11]) ∇u +(Subscript[c, 12])∇v)= (Subscript[g, 1]) - (Subscript[q, 11]) u  - (Subscript[q, 12]) v 
(OverVector[n])·((Subscript[c, 21]) ∇u +(Subscript[c, 22])∇v)= (Subscript[g, 2]) - (Subscript[q, 21]) u  - (Subscript[q, 22]) v

Transient systems have the following addition

(Subscript[m, 11])((∂^2)/(∂(t^2)))u+(Subscript[m, 12])((∂^2)/(∂(t^2)))v+(Subscript[d, 11])(∂/(∂t))u+(Subscript[d, 12])(∂/(∂t))v+∇·...
(Subscript[m, 21])((∂^2)/(∂(t^2)))u+(Subscript[m, 22])((∂^2)/(∂(t^2)))v+(Subscript[d, 21])(∂/(∂t))u+(Subscript[d, 22])(∂/(∂t))v+∇·...

Structural Mechanics

mesh=ToElementMesh[FullRegion[2],{{0,5},{0,1}},"MaxCellMeasure"-> 0.05];
mesh["Wireframe"]

Material Data

PlaneStress[{Y,ν},{u[x,y],v[x,y]},{x,y}]//MatrixForm

youngsModulus=10^3;
poissonRatio=33/100;

op=PlaneStress[{youngsModulus,poissonRatio},{u[x,y],v[x,y]},{x,y}];
pde=op=={0,0};

Multi dependent specification

(Subscript[Γ, D])={DirichletCondition[{u[x,y]==0.,v[x,y]==0.},x==0],DirichletCondition[{u[x,y]==0.,v[x,y]==-1.},x==5]}

This solves the equation.

{uif, vif}=NDSolveValue[{pde,(Subscript[Γ, D])},{u,v},{x,y}∈mesh]

Visualization

Show[{
mesh["Wireframe"[ "MeshElement"->"BoundaryElements"]],
ElementMeshDeformation[mesh,{uif,vif} ]["Wireframe"["ElementMeshDirective"->Directive[EdgeForm[Red],FaceForm[]]]]}]

GraphicsGrid[{{ContourPlot[uif[x,y],{x,y}∈mesh,AspectRatio->Automatic]},{ContourPlot[vif[x,y],{x,y}∈mesh,AspectRatio->Automatic]}},ImageSize->Large]

Structural Mechanics - Boundary Load

mesh["Wireframe"]

op

pde=op=={0,NeumannValue[-1,x==5]};

Multi dependent specification

(Subscript[Γ, D])={DirichletCondition[{u[x,y]==0.,v[x,y]==0.},x==0]}

This solves the equation.

{uif, vif}=NDSolveValue[{pde,(Subscript[Γ, D])},{u,v},{x,y}∈mesh]

Visualization

Show[{
mesh["Wireframe"[ "MeshElement"->"BoundaryElements"]],
ElementMeshDeformation[mesh,{uif,vif} ]["Wireframe"["ElementMeshDirective"->Directive[EdgeForm[Red],FaceForm[]]]]}]

GraphicsGrid[{{ContourPlot[uif[x,y],{x,y}∈mesh,AspectRatio->Automatic]},{ContourPlot[vif[x,y],{x,y}∈mesh,AspectRatio->Automatic]}},ImageSize->Large]

This solves the equation.

{uif, vif}=NDSolveValue[{op=={NeumannValue[-100,x==5],NeumannValue[-1,x>3&&y==1]},(Subscript[Γ, D])},{u,v},{x,y}∈mesh]

Visualization

Show[{
mesh["Wireframe"[ "MeshElement"->"BoundaryElements"]],
ElementMeshDeformation[mesh,{uif,vif} ]["Wireframe"["ElementMeshDirective"->Directive[EdgeForm[Red],FaceForm[]]]]}]

Neumann Values in Systems of Equations

Why are NeumannValues part of the equation?

Consider a system of 2 equations for dependent variables u and v.

((∂^2)/(∂(t^2)))u+(∂/(∂t))u+∇·(-(Subscript[c, 11]) ∇u -(Subscript[c, 12])∇v -(Subscript[α, 11]) u +(Subscript[α, 12]) v+ (Subscript[γ, 1]) ) +(Subscript[β, 11])·∇u +(Subscript[β, 12])·∇v + (Subscript[a, 11]) u+ (Subscript[a, 12]) v == (Subscript[f, 1])
((∂^2)/(∂(t^2)))u+(∂/(∂t))u+∇·(-(Subscript[c, 21]) ∇u -(Subscript[c, 22])∇v -(Subscript[α, 21]) u +(Subscript[α, 22]) v+ (Subscript[γ, 2]) ) +(Subscript[β, 21])·∇u +(Subscript[β, 22])·∇v + (Subscript[a, 21]) u+ (Subscript[a, 22]) v == (Subscript[f, 2])

In this case, the generalized Neumann conditions are

(OverVector[n])·((Subscript[c, 11]) ∇u +(Subscript[c, 12])∇v)= (Subscript[g, 1]) - (Subscript[q, 11]) u  - (Subscript[q, 12]) v 
(OverVector[n])·((Subscript[c, 21]) ∇u +(Subscript[c, 22])∇v)= (Subscript[g, 2]) - (Subscript[q, 21]) u  - (Subscript[q, 22]) v

Mini FEM derivation (single stationary PDE):

(∫_Ω)∇·(-c ∇u  -α u  + γ ) ϕ+β·∇u ϕ+ a u ϕ - f ϕdx=0 .

Integration by parts gives:

(∫_Ω)(c ∇u  +α u  - γ  ) ·∇ϕ+β·∇u ϕ+ a u ϕ - f ϕ dx - (∫_(∂Ω))(OverVector[n])·(c ∇u+α u  - γ )ϕdx=0 .

The generalized Neumann conditions replace the normal derivative (comes from FEM derivation) (OverVector[n])·((Subscript[c, 11]) ∇u +(Subscript[c, 12])∇v) with the value (Subscript[g, 1]) - (Subscript[q, 11]) u  - (Subscript[q, 12]) v

(* plane stress *)
Thread[PlaneStress[{1,1/3},{u[x,y],v[x,y]},{x,y}]=={0,NeumannValue[-1,x==5]}//Activate]//TableForm

If NeumannValue[-1,x==5] floats around in the list of equations, we do not know with which equation to associate it!

Structural Mechanics + Thermo = Thermoelasticity

Thermocouple

Show[
Graphics[{Text["Aluminum",{len 4/5,3hi}],Text["Carbon steel",{len 4/5,-hi}],Text["100C",{-0.01,hi}]}]
,Graphics[Line[{{0,-hi},{0,3 hi}}]],Graphics[Table[Line[{{0,y},{-0.0025,y-0.0025}}],{y,3 hi,-hi,-hi/2}]],mesh["Wireframe"["MeshElementStyle"->{Directive[EdgeForm[Green],FaceForm[Green]],Directive[EdgeForm[Red],FaceForm[Red]]}]]
,mesh["Wireframe"["MeshElement"->"BoundaryElements"]]
]

Geometry

len=0.1;
hi=0.005;
bmesh=ToBoundaryMesh[
"Coordinates"->{{0,0},{len,0},{len,hi},{len,2 hi},{0, 2hi},{0,hi}},"BoundaryElements"->{LineElement[{{1,2},{2,3},{3,4},{4,5},{5,6},{6,1},{3,6}}]}]

bmesh["Wireframe"]

Add region markers for Aluminum and Carbon Steel

mesh=ToElementMesh[bmesh,"RegionMarker"->{{{len/2,3 hi/2},1},{{len/2,hi/2},2}}]

mesh["Wireframe"["MeshElementStyle"->{Directive[FaceForm[Green]],Directive[FaceForm[Red]]}]]

Boundary conditions: Fixed constraints and temperature at the left hand side :

Γ1={DirichletCondition[{u[x,y]==0,v[x,y]==0,Θ[x,y]==100},x==0]};

Boundary conditions: Convective flux over lover and upper edges:

(Subscript[T, ∞])=15; (* Temperature at ∞ in C *)
(Subscript[T, ref])=0; (* Reference temperature in C*)
h=5; (* free boundary convection (W/m^2C) *)
Γ2=NeumannValue[h Θ[x,y]-(Subscript[T, ∞]),y==0||y==2 hi];

PDE model:

pde={
PlaneStress[{YY,νν},{u[x,y],v[x,y]},{x,y}]-te.Inactive[Grad][Θ[x,y],{x,y}]+te*(Subscript[T, ref])=={0,0},
Diffusion[κ,Θ[x,y],{x,y}]==Γ2
};

Thermophysical material specification for each sub domain:

(Subscript[TC, Al])=83;(Subscript[TC, CS])=111; (* Thermal conductivity (W/(m*K)):*)
(Subscript[Y, Al])=70*10^9;(Subscript[Y, CS])=210*10^9; (* Young modulus (Pa) *)
(Subscript[ν, Al])=0.33;(Subscript[ν, CS])=0.29; (* Poisson's ratio *)
(Subscript[TE, Al])=23*10^-6;(Subscript[TE, CS])=12*10^-6;(* Thermal expansion (1/K)*)

PDE coefficient specification:

YY=If[y≤hi,(Subscript[Y, Al]),(Subscript[Y, CS])];
νν=If[y≤hi,(Subscript[ν, CS]),(Subscript[ν, CS])];

κ=If[y≤hi,(Subscript[TC, Al])*IdentityMatrix[2],(Subscript[TC, CS])*IdentityMatrix[2]];

psFactor=Y/(1-ν^2); (* This is for PlaneStress *)
te=If[y≤hi,Evaluate[(psFactor /.{Y->(Subscript[Y, Al]),ν->(Subscript[ν, Al])})* {(Subscript[TE, Al]),(Subscript[TE, Al])}],Evaluate[(psFactor/.{Y->(Subscript[Y, CS]),ν->(Subscript[ν, CS])}) *{(Subscript[TE, CS]),(Subscript[TE, CS])}]];

Solve the equation:

{uif,vif,Θif}=NDSolveValue[{pde,Γ1}, {u,v,Θ},{x,y}∈mesh]

Visualization and post processing:

Show[{
mesh["Wireframe"["MeshElement"->"BoundaryElements"]],
ElementMeshDeformation[mesh,{uif,vif},"ScalingFactor"->5]["Wireframe"["ElementMeshDirective"->Directive[EdgeForm[Red],FaceForm[]]]]}]

ContourPlot[uif[x,y],{x,y}∈mesh]

ContourPlot[vif[x,y],{x,y}∈mesh]

ContourPlot[Θif[x,y],{x,y}∈mesh]

Θif[len,hi]

Interface - API

Typical problem setup:

Ω=ImplicitRegion[0≤x≤5&&0≤y≤1,{x,y}]

Μ=ToElementMesh[Ω,"NodeReordering"->True]

Γ={DirichletCondition[u[x,y]==0, x==0],
DirichletCondition[u[x,y]==1,x==5]}

pde=-Div[Grad[u[x,y],{x,y}], {x,y}]-1/2

uif=NDSolveValue[{pde==0,Γ},u,{x,y}∈Μ]

Plot3D[uif[x,y],{x,y}∈Μ]

Low level access:

{nds}=NDSolve`ProcessEquations[{pde==0,Γ},u,{x,y}∈Μ];

nds

solutionData=nds["SolutionData"][[1]]

nds["FiniteElementData"]["Properties"]

You could use NDSolve`Iterate and NDSolve`ProcessSolution or, at a lower level, extract data via “Methods” or “Properties”

methodData=nds["FiniteElementData"]["FEMMethodData"]

Get the PDE coefficients

pdeData=nds["FiniteElementData"]["PDECoefficientData"]
pdeData["All"]

Get the stationary system matrices

dPDE=DiscretizePDE[pdeData,methodData,solutionData]

(*{l,s,d,m}=dPDE["SystemMatrices"]*)
{l,s,d,m}=dPDE["All"]

MatrixPlot[s]

Get the discretized boundary conditions

bcData=nds["FiniteElementData"]["BoundaryConditionData"]

dBC=DiscretizeBoundaryConditions[bcData,methodData,solutionData]

Deploy the boundary conditions

DeployBoundaryConditions[{l,s,d,m},dBC]

Solve

result=LinearSolve[s,l];

Get the mesh

nr=NDSolve`SolutionDataComponent[solutionData,"Space"]

mesh=nr["ElementMesh"]

Visualize

ElementMeshPlot3D[result,mesh]

Interface - API - Documentation

All low level functions are documented!

Initialization
InitializePDECoefficients — initialize partial differential equation coefficients
InitializeBoundaryConditions — initialize boundary conditions
InitializePDEMethodData — initialize partial differential equation method data

PDECoefficientData ▪ BoundaryConditionData ▪ FEMMethodData

Discretization
DiscretizePDE — discretize initialized partial differential equations
DiscretizeBoundaryConditions — discretize initialized boundary conditions

DiscretizedPDEData ▪ DiscretizedBoundaryConditionData

Solution
DeployBoundaryConditions — deploy discretized boundary conditions into discretized partial differential equations
LinearSolve — solve linear systems of equations

Tutorial

Finite Element Programming

Guide page

Finite Element Method Guide

Thanks and - NDSolve - Rocks!

Summary:

Linear, variable coefficient, transient systems of PDEs over arbitrary regions in
1D/2D/3D

Init

Needs["NDSolve`FEM`"]
$HistoryLength=0;

SetOptions[ContourPlot,{ColorFunction->"TemperatureMap",AspectRatio->Automatic}];
SetOptions[RegionPlot,AspectRatio->Automatic];
SetOptions[#,{ColorFunction->"TemperatureMap",AspectRatio->Automatic,Boxed->False}]&/@{Plot3D,Visualization`Core`Plot3D};

(*ApplyProperty[a_List,property_String]:=#[property]&/@a*)

ClearAll[insertMaterial]
insertMaterial[val_,p_Piecewise]:=Module[{materialRules,conditions,defaultMaterialRule},materialRules=p[[1,All,1]];
conditions=Join[p[[1,All,2]],{True}];
defaultMaterialRule=p[[2]];
Internal`ToPiecewise[conditions,Join[val/.materialRules,{val/.defaultMaterialRule}],False]]


ClearAll[mkRulePiecewise]
mkRulePiecewise[ruleLhs:{__Symbol},a_List]:=MapThread[Rule,{ruleLhs,a}]
mkRulePiecewise[ruleLhs:{__Symbol},p_Piecewise]:=Module[{conditions,rules},conditions=Join[p[[1,All,2]],{True}];
rules=MapThread[Rule,{ruleLhs,#}]&/@Join[p[[1,All,1]],{p[[2]]}];
Internal`ToPiecewise[conditions,rules,False]]

ClearAll[mkMaterialDispatcher]
mkMaterialDispatcher[ruleLhs:{__Symbol},input_]:=Block[{newPW},newPW=mkRulePiecewise[ruleLhs,input];
If[Head[newPW]===List,With[{newPW=newPW},Function[in,in/.newPW]],With[{newPW=newPW},Function[in,insertMaterial[in,newPW]]]]]

Diffusion[k_,u_,X:{x__Symbol}]/;!ArrayQ[u]:=With[{res=Catch[iDiffusion[k,u,X,Diffusion],NDSolve`FEM`FEMError[___]]}, res/;Unevaluated[res]=!=$Failed
]

ClearAll[iDiffusion];
iDiffusion[kIn_,u_,X:{x__Symbol},msghead_Symbol]/;!ArrayQ[u]:=Module[{k,a},
k=kIn;
If[NumberQ[kIn],k=kIn*IdentityMatrix[Length[X]]];
 If[VectorQ[kIn,NumericQ],
NDSolve`FEM`FailIf[Length[kIn]=!=Length[X], msghead,"Input vector `1` needs to be of the same length as the nuber of independent variables (`2`)",kIn,Length[X]];
k=DiagonalMatrix[kIn];
];
a=PiecewiseExpand[Piecewise[{{k,True}}]];
Inactive[Div][a.Inactive[Grad][u,X],X] 
]

PlaneStress[{Y_,nu_},U:{u_,v_},X:{x_,y_}]:=With[{res=Catch[iPlaneStress[PiecewiseExpand[Piecewise[{{{Y,nu},True}}]],U,X,PlaneStress],NDSolve`FEM`FEMError[___]]},res/;Unevaluated[res]=!=$Failed]

ClearAll[iPlaneStress];
iPlaneStress[a_,{u_,v_},X:{x_,y_},msghead_Symbol]/;(Head[a]==List&&Length[a]==2)||(Head[a]==Piecewise):=Module[{pStress,Y,nu,valF},pStress=-Y/(1-nu^2)*{{{{1,0},{0,(1-nu)/2}},{{0,nu},{(1-nu)/2,0}}},{{{0,(1-nu)/2},{nu,0}},{{(1-nu)/2,0},{0,1}}}};
If[Head[a]==List&&Length[a]==2&&NumberQ[a[[2]]],NDSolve`FEM`FailIf[TrueQ[Sqrt[a[[2]]]==1],msghead,"`1` must not be `2`","nu^2","1"];];
valF=mkMaterialDispatcher[{Y,nu},a];
{Inactive[Div][valF[pStress[[1,1]]].Inactive[Grad][u,X],X]+Inactive[Div][valF[pStress[[1,2]]].Inactive[Grad][v,X],X],Inactive[Div][valF[pStress[[2,1]]].Inactive[Grad][u,X],X]+Inactive[Div][valF[pStress[[2,2]]].Inactive[Grad][v,X],X]}]

(*
PlaneStrain[{Y_,nu_},U:{u_,v_},X:{x_,y_}]:=With[{res=Catch[iPlaneStrain[PiecewiseExpand[Piecewise[{{{Y,nu},True}}]],U,X,PlaneStrain],NDSolve`FEM`FEMError[___]]},res/;Unevaluated[res]=!=$Failed]

ClearAll[iPlaneStrain];
iPlaneStrain[a_,{u_,v_},X:{x_,y_},msghead_Symbol]/;(Head[a]==List&&Length[a]==2)||Head[a]==Piecewise:=Module[{pStrain,Y,nu,valF},pStrain=-((1-nu)*Y)/((1-2*nu)*(1+nu))*{{{{1,0},{0,(1-2*nu)/(2*(1-nu))}},{{0,nu/(1-nu)},{(1-2*nu)/(2*(1-nu)),0}}},{{{0,(1-2*nu)/(2*(1-nu))},{nu/(1-nu),0}},{{(1-2*nu)/(2*(1-nu)),0},{0,1}}}};
If[Head[a]==List&&Length[a]==2&&NumberQ[a[[2]]],NDSolve`FEM`FailIf[TrueQ[Abs[a[[2]]+1]≤$MachineEpsilon||Abs[a[[2]]-1/2]≤$MachineEpsilon],msghead,"`1` must not be `2`","nu","-1 or 1/2"];];
valF=mkMaterialDispatcher[{Y,nu},a];
{Inactive[Div][valF[pStrain[[1,1]]].Inactive[Grad][u,X],X]+Inactive[Div][valF[pStrain[[1,2]]].Inactive[Grad][v,X],X],Inactive[Div][valF[pStrain[[2,1]]].Inactive[Grad][u,X],X]+Inactive[Div][valF[pStrain[[2,2]]].Inactive[Grad][v,X],X]}]*)